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We present two methods for approximating the mapping 
between two systems exhibiting generalized synchronization. 
If the equations of motion are known then an analytic approx- 
imation to the mapping can be found. If time series data is 
used then a numerical approximation can be found. 
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The subject of synchronization between identical sys- 
tems (denoted here by IS) has been of interest since the 
time of Huygens. Over the last decade it has become 
clear that even chaotic systems can be synchronized [[lj. 
One example is drive-response synchronization, where 



~dl 
dy 
dt 



F(x) 

F(y)+E(x,y). 



Here, E(x, y) denotes coupling between the drive system 
(x) and the response system (y). If F is deterministic, 
and if E(x,x) = 0, then the systems are synchronized 
if y(i*) = x(i„). Because of determinism this condition 
remains true for t > t*. 

Recently, papers discussing a more general idea of 
synchronization have appeared in the literature. Drive- 
response dynamics for this type of synchronization is 
given by, 



rfx 
~dt 
dy 
dt 



F(x) 
G(y;x) 



(1) 



where G and F are permitted to be different functions. 
In principle, x £ M d , y £ M r , and the dynamics takes 
place in M d+r . Intuitively, Generalized Synchronization 
(GS) occurs if the response state, y, is related to the drive 
state, x, by a time independent function, y = 0(x). If GS 
occurs then the dynamics takes place ona<i dimensional 
invariant manifold in lR d+r . 

Much of the work on GS has focused on three areas. 
The first area focuses on defining GS. Various definitions 
have been proposed ||-||] . Reference Q suggest that sub- 
tleties associated with unstable periodic orbits imply that 
more that one definition may be required. The second 
area focuses on mathematical properties of </>. Rigorous 
results about the smoothness of </>, and the relationship 
between smoothness and Lyapunov exponents exist ||^| . 
Also, numerical methods for determining the properties 



of <fi exist pj. Since GS has been observed in experimen- 
tal systems it is structurally stable. Mathematical lit- 
erature regarding the existence, stability and smoothness 
of invariant manifolds is also relevant 0. The last major 
area of research has focused on detecting GS from time 
series data §|uj. The methods are indirect in the sense 
that they cither do not approximate <p, or the approxima- 
tions are local (often resulting in as many approximations 
as date points). 

This manuscript opens a new research direction. In- 
stead of seeking properties of 0, or indirect evidence of 
its existence, we believe it is better to go after the func- 
tion itself. Therefore, we present methods for analytically 
and/or numerically constructing a single smooth function 
which globally approximates d>. If the equations of mo- 
tion are known then an analytic approximation for <f> can 
be obtained. (To our knowledge, this is the only tech- 
nique for analytically approximating </>.) The numerical 
method uses time series from the two systems to calcu- 
late a statistic which can be used to infer the existence 
of stable GS. The numerical method also gives a global 
approximation for the function, y = </>(x). (We argue, 
implicitly, that if <fi and/or </> -1 exist but are not well ap- 
proximated by smooth functions then their usefulness is 
limited since their mathematical properties are probably 
"so bad" they prohibit most applications of GS.) 

An important application for GS comes from control 
theory. Typically, control schemes work better when the 
complete state of the plant is known. The application 
uses measurements from the plant (F) as drive input to 
an approximate model of the plant (G). If GS occurs 
then the state of the plant can be approximated from 
the state of the model via x = 4>~ l (y). This, and most 
other applications, require a stable GS manifold. 

Recently, several criteria have appeared for designing 
coupling which results in a stable IS manifold ]TT] , p^[ . 
We report here that it is straightforward to show that a 
criteria for linearly stable GS is jl3| 

A^(D y G[0(x);x]) 
-R[A X ] > (HP" 1 [D y G[0(x);x] - (D y G)] P||) . 

Here, (•) denotes a time average over the driving trajec- 
tory, 5i[Ai] is the eigenvalue of A with the largest real 
part, and P is a matrix whose columns are the eigenvec- 
tors of A. Also, D y G denotes the Jacobian of G with 
respect to y. This criteria implies that if 4> an d 4>~ l are 
known then one can estimate the state of the plant (x) 
from the state of the model (y) by design coupling which 
guarantees stable GS. 



The analytical method used to approximate <fi is based 
on approximating center manifolds. Although the appli- 
cation to GS is new, complete discussions about approx- 
imating center manifolds (with examples) can be found 
in many text books JT3]. Therefore, our discussion will 
be brief. 

Assume the drive and response systems are given by 
Eq. (Q). Taking the total time derivative of y = 0(x), 
and using Eq. (|lj), implies that 



G[0(x);x]-[D x 0(x)]-F(x)=O 



(2) 



on the synchronization manifold. Here D x </> is the Jaco- 
bian of (p. Equation (J2j) is interpreted as a partial dif- 
ferential equation for the unknown function, 0(x). The 
same t ype of equation arises when estimating center man- 
ifolds g. 

Typically, Eq. (g) can not be solved exactly. Therefore, 
approximate the solution by the series </>(x) = A + B • 
x + x- M- x+ .... Next, insert the series into Eq. (g) 
and rewrite the results as a polynomial in powers of x. 
The coefficients of this polynomial are functions of the 
parameters of F and G as well as the elements of A, B, 
M, etc. Also, this polynomial must hold for all x on the 
driving trajectory. If this trajectory is not a fixed point 
then it is reasonable to assume that the polynomial can 
hold only if the coefficient of each power of x vanishes. 
By equating each coefficient to zero we form a set of 
algebraic equation involving the parameters of F, G, and 
the elements of A, B, M, etc. The approximation to 
0(x) is obtained by solving these algebraic equations for 
A, B, M, etc in terms of the parameters of F and G. 

Although conceptually straightforward, performing 
this procedure on anything but the simplest examples is 
very tedious, and soon grows beyond what can be done 
by hand. However, these calculations are not beyond the 
power of modern symbolic manipulation software. In- 
deed, the results presented below were obtained using 
MAPLE . It was relatively straightforward to write a 
MAPLE program which produced these answers. Once 
the program was written the total run time was less than 
10 minutes. 

The approximation that one obtains for the GS mani- 
fold should hold near the attractor for the drive dynam- 
ics, however, it is not likely to be globally well defined. 
Although the results will not be presented, we have used 
a similar analysis to approximate x = </> _1 (y) for all of 
the examples discussed below. 

If the GS manifold is stable then we can numeri- 
cally approximate <fi from time series data. The nu- 
merical method used to approximate <fi is similar to one 
used by several authors to make empirical global mod- 
els from time series data. Begin by assuming one has 
two data sets, x(rtAi) e R d and y(nAi) € R r , with 
n = 1,2, ... ,N, which represent simultaneous measure- 
ments of the drive and response systems at a sampling 
rate At. (If necessary, vector representations of the dy- 
namics can be obtained from scaler time series via em- 



bedding techniques |L(|.) A measure for the dynamics of 
the drive system can be approximated by [fL6| 



1 



JV 



p(z) = hm — > 8\z — x(n)l 

n=l 

Since the exact functional form of <f> is unknown the 
best one can hope for is a series expansion 



K 



4>(z) = lim > p 



(z). 



(3) 



1=0 



Here, the pW's are r dimensional expansion coefficients, 
which must be determined, and the 7i"W(z)'s represent 
some set of basis functions. Several authors have demon- 
strated the advantage of using a basis set which is or- 
thonormal on p(z) , and they show how to construct such 
a basis from data using Gramm-Schmidt @@. The 
summation index, I is used to identify the individual ba- 
sis functions. 

Once the basis set has been constructed, each expan- 
sion coefficient, p*- 1 * 1 , can be obtained by multiplying both 
sides of Eq. (g|) by 7r®(z) p(z) and integrating over all 
space. Because of the orthonormality of the basis set we 
obtain 



lim 

iV-> oo N 



1 N 



r(I) 



x(n)], 



(4) 



n=l 



where we have used y(n) = 0[x(n)] on the GS manifold. 
Thus, Eqs. (||) and (|j) are used to approximate <fi from 
time series x(n) and y{n). 

The last task is to determine the order at which to 
truncate the series in Eq. (||) so as to not over fit the data. 
This is done by using the minimum description length 
(MDL) criteria. This criteria is similar to the maximum 
likelihood principle associated with least squares fitting 
of data |fiq , ^9| . However, unlike maximum likelihood, 
MDL is capable of determining the optimal order at with 
to truncate Eq. (g). The MDL function we use is given 
by 

rN 



Xmdl = — [In (27ra 2 ) + l] 



+ N V 



Mi) 



K 



1=0/3=1 



(See Ref. |l1| for a complete derivation of this function.) 
Except for a positive constant (which we neglect), the 
first term is the usual prediction error from the maxi- 
mum likelihood principle. Indeed, a 2 is the least squares 
prediction error obtained when predicting the y(n)'s from 
the x(n)'s. 

The remaining terms are penalties which increase as 
more terms in Eq. (g) are retained and the model be- 
comes more complex. N p is the total number of nonzero 



p(i)' s retained in Eq. (|^). In our implementation, a com- 
ponent of pW's is set to zero if its statistical significance 
is not distinguishable from zero 8f is the relative 

accuracy of the (3 component of p^, rj is the relative 
accuracy of a 2 , and 7 = 32. 

To illustrate the analytical and numerical techniques 
we applied them to examples using the Lorenz equations 



— = s(x 2 - Xl ) 



dt 



rxi - xix 3 - x 2 



(5) 



dx 3 

— — = x x x 2 - bx 3 , 
dt 



as the drive system. 

The coupling between drive and response systems usu- 
ally involves one of two cases. The first case arises when 
the physical processes responsible for the coupling are 
known so one has an explicit equation for the coupling. 
For this case one solves Eq. (||) as discussed above. Below 
we consider the second case where the response system is 
given by G(y) 4- E[0(x) — y]. Here, the coupling obeys 
E(0) = and is used to insure that the GS manifold is 
stable. The problem with this case is that we can not 
evaluate E[0(x) — y] because we do not know, a priori, 
the form of </>(x). For the examples discussed below this 
problem is overcome by calculating <p in two stages. 

In the first stage we calculate <p using diagonal cou- 
pling E[/0(x) — y] = e['0(x) — y] where i/> is an arbitrarily 
function. The <p calculated in this first stage clearly de- 
pends on tp. In the second stage we force tp(x.) = 0(x). 
This second stage insures that E[i/>(x) — y] = on the 
GS manifold. 

Two trivial tests of the analytic method involved defin- 
ing y = </>(x) = [x\ +ax2 + (3x2, x 2, x 3 ] for one test, and 
y = 0(x) = [x\ + ax 2 , [3x2, x 3 + ix 2 ] f° r the other. For 
each test we obtain a response system, y = G(y), by 
taking the time derivative of y, using Eq. (||) to resolve 
the vector field, G(y), and adding diagonal coupling. For 
these test, the response systems are the Lorenz system 
after a nonlinear change of coordinates, and the analytic 
method easily recovered the GS manifolds, y = </>(x). 

A final test of the analytic procedure used the following 
response system 



— = r(l + A) yi - Vl y 3 - y 2 
^ =ViV2 -b(l + r))y 3 , 



(6) 



with diagonal as discussed above. 

For this example we could only approximate <p. The 
approximation contained three arbitrary constant, thus it 
is not unique. We selected values for two of them so that 
<p has a simple form. (For the trivial examples discussed 



above this choice always lead to the the "correct" equa- 
tion for </>(x)). The third constant appears trivially in 
B 33 , is of order (A, 5, 77), and is denoted by K below. Fi- 
nally, the approximation is simplified by retaining terms 
that are second order in x, first order in (5, A, rf), and in 
the limit of large coupling strength, e. (Thus, we exam- 
ine a case where the response system is close to the drive 
system.) 

With these criteria in mind we found that </3(x) is given 
by A = 0, 



B-l + i 



-(r 2 A-s 2 <5) -s{r+l-s)S 
-(r(s - 1)A - s 2 S) (r 2 A- 




>5) 



where 1 is the identity matrix, and the three tensor M 
is given by M (1) = M (2) = 0, and 



M 



(3) 



1 

r 



r(s-l)A-s 2 S 
b-2s 








s(b-2s)(r+l-s)S q 




fc-2 





In these equations, T = (2r 2 + 3s 2 — 2,s + l). Furthermore, 
it clear that this transformation satisfies <fi = 1 in the 
limit 5, A, 77 — > 0. 

To test the numerical method we first demonstrate 
that it can determine the correct form of (p for stable 
GS from time series data. To accomplish this we used 
Eq. (||) (with ,s = 16, b = 4, and r = 46) as the drive 
system and a response system obtained from y = </3(x) = 
[xi — 0.01a; 2 ;, 0.95x2, x 3 + 0.03x 2 ]. The systems were cou- 
pled via the 7/2 equation using e[(0.95i£2+noise)— y 2 ]- The 
noise was Gaussian white with zero mean and standard 
deviation, 15a. Here, 15 is approximately the standard 
deviation of X2, and a = or 0.05. e = 10 was used be- 
cause, with ?/2 coupling and a chaotic driving trajectory, 
IS is "stable" for e > 4 0. 

The numerical procedure was given N — 4000 simulta- 
neously recorded values of x and y at a sampling interval 
of At = 0.02. The results (see table |) indicate that the 
numerical procedure found a good approximation to (p, 
even in the presence of small amounts of noise. 

To further test the numerical method we used Eqs. (||) 
and (§) (the same values for s, b, and r) and a drive 
signal, e[(x2 + noise) — 2/2], coupled to the 7/2 equation. 
These tests used simultaneously recorded scalar time se- 
ries of the same length and sampling interval given above. 
Scalars were obtained using the arbitrarily chosen projec- 
tions 

Sd(n) = xi(n) — 2.5x2(n) + 0.75x 3 (n) 
s r {n) = ~0.byi{n) + 1.5y 2 {n) - y 3 {n). 

Each scalar time series was independently rescaled to 
mean zero and standard deviation one, and an attrac- 
tor for each time series was reconstructed using a time 
delay embedding |lq|. 



Factor 


01 


(7 = 

02 


03 


<h 


■j = 0.05 

02 


03 


const 


0.00101 


0.00055 


-0.00194 


0.0822 


0.00215 


0.00475 


X\ 


1.00 








1.023 





-0.0390 


x 2 





0.950 





-0.0107 


0.954 


0.0147 


X3 








1.00 


-0.0061 





1.003 


XiXi 











0.000982 





0.000166 


X\X2 











-0.000304 





-0.000259 


X1X3 











-0.000233 





0.000381 


x\ 








0.0300 








0.0297 


X2X3 




















xl 


-0.0100 








-0.00994 









TABLE I. Numerical approximations for the transforma- 
tion 0i = xi — 0.01x|, <j)2 = 0.95:E2, and c/>3 = xz + 0.03a;|. If 
the calculate value of (f>j was of order 10 -5 or less then it was 
set to zero. 
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FIG. 1. The sudden drop at e ~ 4 indicates the onset of sta- 
ble synchronization and stable generalized synchronization. 



The results of our attempts to approximate <fi for S = 
A = 77 = (IS) and S = 0.02, A = 0.04, rj = -0.03 (GS) 
are shown in Fig. ||. The figure shows that Xmdl expe- 
riences a sharp drop at e ~ 4 when the drive/response 
systems are identical and a less sharp drop for GS. The 
drop implies that the numerical procedure has found a 
relatively accurate approximation for y = c/>(x), so the 
GS manifold is stable. Also, the figures shows that the 
procedure deteriorates gracefully in the presence of noise. 
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In conclusion, we have presented an analytical and a 
numerical method for approximating the mapping that 
defines the invariant manifold associated with generalized 
synchronization. The author would like to thank Drs. N. 
F. Rulkov, L. M. Pecora, B. R. Hunt, and J. Stark for 
valuable discussions and comments that lead to this work. 



